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Abstract 

Theory of solutions in energy representation (ER method) developed by Matubayasi 
and Nakahara provides with an approximate way of calculating solvation free energies 
(or, identically, the excess chemical potentials) from atomistic simulations. In this doc¬ 
ument we provide some derivation details of this, to our opinion, theoretically involved 
method, which will help a non-specialist to follow. There are three points which differ 
this document from a regular textbook on statistical mechanics or research articles: 

1. Derivation is detailed and all approximations are explicitly stated; 

2. Statistical mechanics derivations are performed in NPT-ensemble; 

3. We perform the derivations for the case when a molecule is represented as a set 
of (atomic) sites interacting via spherically symmetric potentials (a classical Force 
Field representation). 

In ER method, a new collective coordinate is introduced - the interaction energy of a 
solute and a solvent molecule. The excess chemical potential is expressed as a functional of 
the solute-solvent density distribution dehned over the collective variable. The functional 
can be approximated by the Percus’s method of functional expansion, which leads to the 
end-point (not dependent on the A-coupling path) free energy expression. 

As a side result, we prove that the solvation free energy is always equivalent to the 
excess (over ideal) chemical potential, and not only at inhnite dilution or when internal 
molecule degrees of freedom are not affected by solvation as it is sometimes wrongly 
believed. 
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1 Introduction 

We provide detailed derivation of the theory of solutions in energy representation (ER method) 
developed by Matubayasi and Nakahara [l|, l2|, y, y, l5[ . ER method provides with an approximate 
way of calculating solvation free energies (or, identically, the excess chemical potentials) from 
atomistic simulations. The method can be seen as bridge between the molecular simulations and 
the classical density functional theory (DFT). It is quite common nowadays to model molecular 
interaction on the level of the classical force field approximation, which implies that a molecule 
is represented as a set of (atomic) sites interacting via spherically symmetric potentials. 

In the first part of the manuscript we define and derive the expression for solvation free 
energy (SEE) in NPT-ensemble for the case of classical force field representation of molecular 
interactions. We prove that SEE is identical to the excess chemical potential. Also, we obtain 
the Kirkwood’s charging formula expressing the excess chemical potential via the solute-solvent 
density distribution. 

In the second part of the manuscript we provide details on some important relations of ER 
method for the case of NPT-ensemble. In ER method, a new collective coordinate is introduced - 
the interaction energy of a solute and a solvent molecule. We show that the Kirkwood’s charging 
formula is valid also for the solute-solvent distribution function in energy representation. Later 
this expression is reformulated as a functional of density distribution. The Percuss method of 
functional expansion is used to obtain the hypernetted chain (HNC)- and Percus-Yevick (PY) 
-like approximations of this functional. The final formula for the excess chemical potential 
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heuristically combines expressions from different approximations and employs different input 
functions. 


2 Excess chemical potential in NPT-ensemble 

2.1 Some definitions 

We consider a system with Ng = 1 solute and N^j solvent molecules in isothermo-isobaric 
ensemble (NPT-ensemble). 

We describe the interactions between the molecules in the force held approximation at 
the level of classical mechanics. Each molecule is represented as a set of atoms (better to 
say interaction sites), which interact with each other via bonded and nonbonded potentials 
present in the given force held (e.g. OPTS, CHARMM, AMBER, etc). Each interaction site 
is considered as a separate object, which has its own translational degrees of freedom and 
translational partition function. Therefore, when we talk about a set of coordinates which 
dehne the position of a molecule Xj we mean the positions of all atoms which belong to this 
molecule, where index a runs over all atoms (nt) of the molecule of type t: 

Xi = {ri,a}"Ll (1) 

where t is the molecule type, e.g. s denotes solute, w denotes solvent, and the coordinates of 
atom a of molecule: 

^i,a Vi,a) ^i,aj 

Each atom has its own momentum; pi^a 

2.2 Parametrized Hamiltonian 

Here and after we mostly adopt the notations used in the Appendix of the Shirts et ah publi¬ 
cation j^. 

The excess chemical potential can be calculated in the process of gradual switching on 
the intermolecular interactions between a solute molecule and the solvent. We introduce A 
parameter which controls the degree of coupling between the solute and solvent molecules, 
such that when A=0 the interactions are absent and when A=1 interactions are at full coupling. 
Since, only solute-solvent interaction potential Usw^x(xs, x^^*) depends on A, the potential energy 
function of the system can be written as follows: 

= d'(x^) + ^Ms^,A(Xs,X^,i) -F 17^^(x^“) (2) 

i=l 

where subscript s denotes solute, subscript w denotes solvent, ^'(xs) is the potential energy of 
the solute molecule, x^, j is the position of solvent molecule, is the number of solvent 
molecules, Usw,x is the A-dependent solute-solvent interaction potential, Uww is the potential 
energy of the solvent molecules, x^*" is the short notation of positions of all solvent molecules. 
The total Hamiltonian can be written as: 

Hx = K{pg,p^'^)+ Ux{xg,x^^) (3) 
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where the kinetic energy is written as: 


where rUs^a and 

ingly. 


K{Ps,P^-) = Y, 


P. 


s,a 


N't!) PI'W 2 
r 


2m. 


EE 


w,t,a 


2m,. 


(4) 


a=l i=l 0=1 

are the masses of atoms of solute and solvent molecules, correspond- 


2.3 Partition functions with non-parameterized Hamiltonian 

Keeping in mind that we consider the system with a single solute molecule Ng = 1, we will 
write explicitly the terms with Ng in the derivations. Later this will help us to show that the 
SFE is always equal to the excess chemical potential. 

2.3.1 Case of solution 

The partition function in NPT ensemble can be written as: 

A{Ng, iV^, P, T) = d e-^^^Q{Ng, K, T) (5) 

where V is an arbitrary constant which makes the partition function dimensionless, f3 is 
{kBT)-\ ks is the Boltzmann constant, Q{Ng, N.,,j,V,T) is the canonical partition function, 
which has the following form: 


Q{Ng,N,,,V,T) = 


“+00 


h^N,n nu, M ! flSNsTls jY I 


dp^-dp 


'V 


<ixf<ixyexp[-/3i/(pf,py 






)] 


where h is the Planck’s constant. Multiplication by h~^ serves as a quantum correction for 
purely classical partition function [^. The factorials of number of atoms in the system appear 
due to indistinguishably of atoms belonging to the molecules of the same type. Each integration 
symbol denotes integration over multiple coordinates. Differential is the short notation 

for Symbol V at the integration sign - fy - reflects that integration limits are 

bound by the system’s volume. 

In the case of classical statistical mechanics the momenta degrees of freedom are independent 
and can be analytically integrated j^: 


“+00 


h 


-3 


dp 




x,a'- 


h‘^ 


- 1.5 


27rmx,akBT 


A 


-3 

x,a 


where Ax,a is the thermal de Broglie wavelength for atom a in molecule of type x. 
Therefore we get: 


\—3Ns \—3 Nil 

1 r 


e(JV„ I/, r) = J] n 


a=l 


a=l 


N„ 


■Z{Ng,N,p,V,T) 


where Z is the conhguration integral of the system: 


Z{Ng,N,p,V,T)= / dxf«dxf-exp[-^[/( 
Jv 




( 6 ) 

(7) 
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The Gibbs free energy is: 


G{Ns,N^,P,T) = -kBT\nA{N,,N^,P,T) 

The chemical potential of solute in the system can be written as: 

A{N,,N^,P,T) 


/i = G{Ns, N^, P, T) - G{N, -GN^,P,T) = -kBT In 




( 8 ) 


2.3.2 Case of ideal gas 

For later derivations we will use the expression for the chemical potential of non-interacting 
solute molecules at given T, V and Ng. Therefore, we derive it here. Firstly, let us hnd the 
conhguration integral for a single solute molecule: 

Z{Ng = = 0,V,T) = f dxsexp[-P^s{xs)] (9) 

Jv 

The potential energy in the force held (FF) representation is a function only of distances 
between particles and does not depend on their absolute positions. Additionally, we consider 
homogeneous liquid phase. These two facts allow us to change the coordinates of the system 
such that one atom of the solute molecule is located in the origin . Coordinates of all solute’s 
atoms are written as (see Eq. [1]): 




Therefore, we may rewrite Eq. O 



dxg exp [-/3T^(xJ] 


'V 


dr's,idr'g^2-dr'g,ns 


s,2) 



= V 


'V 


dr's^2-dr's,ns^^P [-^^.( 0 , 


-s,2’ •••? ^ 


s,ns 


= V drs, 2 ...drs,n, exp [-/3Ts(0,r5,2, ■•■,rs,nj] = 

Jv 

where we, hrstly, integrated out the position of the hrst atom of the solute molecule which 
released the volume, and, secondly, we dropped the / marks for simplicity. 

To proceed we introduce the following approximation. The bonded potentials in FF rep¬ 
resentation do not allow atoms belonging to the same molecule to move far from each other. 
Therefore, the limits of integration for the rest of solute’s atoms with very high accuracy can 
be reduced to a small volume around the hrst atom (we denote it as Vg). Please note that, 
since we consider a single molecule here these do not affect the combinatorial prefactor of the 
canonical partition function. Therefore we write: 


dx^exp [-/3T^(xJ] = V 


'V 


'Vs 


'Vs 


drg^2---drg^ns exp [-/3T^(0, G,nJ] 


For simplicity we will use the following notations: x* = {r^ 2 , correspondingly, 

dx* = drg^ 2 ---drg^ns- With these notations we have: 



dx^exp [-/3T^(xJ] = V 



dx* exp [-/3T,(0, x*)] = V ■ g,(T) 


( 10 ) 
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where we introduced new function qs{T), which is in some sense corresponds to the internal 
partition function of a solute molecule in FF representation: 


Qs{T) = f dx*exp [-/3d',(0,x*)] (11) 

However, one should note that in this dehnition all conformations of the molecule are taken 
into account in contrast to the usual q{T) dehnition based on vibrational, rotational, electronic, 
etc. partition functions, which are dehned for a single molecular conformation j^. 

With the help of Eq. [10] we may write the conhguration integral (Eq. [7|) for non-interacting 
solute molecules as: 


Z,a{N,,N^ = 0,V,T)= I dxf^exp 


>v 


Ns 




i=l 


'V 


Ns 


<ix.expM>I>(x,)] ) 

( 12 ) 


And the corresponding canonical partition function (see Eq. [6|) is: 

ris j^- 3 Ns 


QuiN,, N„, = 0. V, T) = J] . ,f-(r) 


a=l 


(13) 


The chemical potential for the ideal gas case at given V is written as: 


, Qu(n.,n^ = o,v,t) , rt;Li^w«--jf-(r) 

l^id = -kBTln ^ =-/cgTIn 




.-3(JVs-l) 

mi. %fciir ■■ 9* (r) 


{Ns-l) 


= —ksTln 


VqsiT) 

Nris 
^ s 



a=l 


(14) 


2.4 Solvation free energy and excess chemical potential 

The solvation free energy (SFE) can be dehned as a reversible work required to switch on the 
interactions between a solute molecule and the rest j^. In NPT ensemble this can be written 
as: 


^Gsoiv = -ksTlxi 


A(iV„iV^,F,r,A = l) 

A(W,iV«i,P,r,A = 0) 


(15) 


where the A in the brackets indicate that the partition functions are written with the A- 
parameterized Hamiltonian (Eq. [3|). 

In the next transformation of Eq. [T5] the factor in the nominator appears because 
of the dissemination process [^. When we write the parameterized Hamiltonian we scale 
the solute-solvent interactions only for one solute molecule, which makes this solute molecule 
distinguishable from the rest. This switch from assimilated and disseminated solute molecule 
changes the combinatorial prefactor of the canonical partition function (Eq. [6]). Inside a single 
molecule we consider atoms being physically different, however atoms of the same type from 
identical molecules are physically identical. Therefore the A”'’ factor appears: 


AG.„/„ = -/crTIu 


a;-a(a„a^,p,t) 


/“ d (f) = 1,N^ = 0, H, T)Q{Ns - 1, A^, H, T) 


(16) 
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where A{Ns, P,T) is the partition function with non-parameterized Hamiltonian and the 
canonical partition function in the denominator factorizes into canonical partition function for 
the single solute molecule and the system with — 1 solvent molecules. 

The denominator in Eq. Hn] can be further simplihed: 

d = 1, iV^ = 0, H, T)Q{N, - 1, H, T) = 

(please, again, note that one solute molecule is not identical to the rest and the corresponding 
combinatorial factor reduces by one:) 



n 


A 


-^Ns 


n 


A 




LJ- (Ns - 1)! 

=1 ^ ® ' 0=1 






0,H,T)Z(iV,-l,iV^,H,T) = 


(17) 

(for a large number of molecules in the system only integration of large system volumes 
contribute to the NPT partition function (Eq. [5]), and therefore integration over small volumes 
comparable to 14 can be safely neglected. Thus, the integrals over small volumes 14 and, 
consequently, Z{Ns = = 0,17, T) become independent on total volume 17. Therefore, 

using Eqs. dUlwe can rewrite Eq. [T7I as:l 


= 




a=l 


A—3(Afs —1) n-u, A —3Vu, ^oo / T/\ /" \ 

n (fvT^ n ^^y I .v.z(N,-i,N„v ,Dj ^ 

(18) 

(multiplication and division of Eq. [TH]by A{Ns — 1, N^, P, T) leads to the following:) 

= ■A(iV,-l,iV^,P,r) 

where 17* is the average volume of the {Ng — 1, N^, P, T) system. 

Finally, we write SEE (Eq. dS]) as: 


AGsoiv — —ksTln 
(using Eq. [8] we get:) 


A{Ng,N^,P,T) 


n: 


[AiNg-l,N^,P,T) qs{T)-V*-UyA- 


-3 

a 


= /i + ksTln 


9 .(r)w.n:'.iA-j 


JVns 


we add fc^Tln where Vi is the mean volume of {Ng, N^) system:) 


qs{T) • El • n ^ii 

m 


= /i + ksTln 


V* 

+ ksT In— 

n 


(19) 


where the last term is the work required for one ideal gas particle to expand the volume from 
17* to 17i. In thermodynamic limit the ratio of two volumes tends to 1 and therefore the term 
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vanishes. The hrst term is the chemical potential of the ideal gas of solute molecules (see 
Eq. [TT|) with the average volume of {Ng, P, T) system Vi. Therefore, Eq. [19] is rewritten as: 

AGsoiv = h - h"'" (20) 

Eq. 120] shows that the SEE is always equal to the excess (over ideal) chemical potential. This 
also proves that the excess chemical potential is the reversible work of switching the solute - 
solvent interactions, as dehned in Eq. [16] Therefore, we may write: 

solv — ^ex (21) 

We would like to note two things here. Firstly, SEE equals to the excess chemical potential 
not only at inhnite dillution, as it is sometimes wrongly believed. 

Secondly, there is a wrong statement (at least for the present case of the FF models of 
molecules) in the book of Ben-Naim (Ref. [^, page 200) that ”... only when is unaffected 
by the solvation process, yUex becomes identical with the solvation Gibbs energy ...”. In our 
derivation of Eq. |2I] we explicitly considered the case when internal degrees of freedom of 
molecules, represented by sets of interaction sites, are coupled to other degrees of freedom. 


2.5 Kirkwood charging formula 

In order to make the forthcoming derivations simpler, from now we explicitly consider the case 
of inhnite diluted ssoluteion: Ng = 1. The point here is that we will express the excess chemical 
potential of solute via the particle density distributions. Therefore, considering many solute 
molecules in the system will require to introduce the solute-solute density distribution, which 
will unnecessarily complicate the derivations. Note, that all the forthcoming derivation can be 
straightforwardly extended to the case of multicomponent solvent (see Ref. (inj|). 

Let us dehne the excess chemical potential for the system with parameterized Hamiltonian 
(Eq. Ej) at a certain A value. With Eqs. ITBl and I2T] we get: 


= -ksTln- 


AiNg,N^,P,T,X) 


A{Ng,N^,P,T,X = 0) 

Since the denominator does not depend on A one can write: 

Jo°°^(w) exp [-/3Uxixg,x^-)] 


dl-l'exp 


= -ksT 


dX " A{Ng,N^,P,T,X) 

With the explicit form of the potential function (Eq. |2|) we have: 


dUx 

dX 


dp, 


etc,A 


dX 


2 = 1 


dug 


',A 7 ^22?,2 


^+oo 


dX 


dx'^dxl^^!^^AA^^£l]^ ^ A(x, - x'JA(x^,i - x'J 


2=1 


where (■)^ denotes ensemble average in isothermo-isobaric condition at given A. We can change 
order of integration and take out the derivative from the ensemble average: 


dp, 


etc,A 


^H-oo 


dX 


dx[dx[, 


dug^^x{^'^,^'J 

dX 






( 22 ) 


In the right hand side there is the pair solute-solvent density distribution in NPT-ensemble 
by dehnition (see e.g. Eq. (2.5.13) of Ref. jsl, but mind that for density distributions of 
non-identical particles the sum should include terms with i = j): 













= / (ix'dx! 




psw,x{^si Xui) 


dX J " " dX 

Finally, the excess chemical potential can be written as an integral over lambda: 




'■+°° du 

dxc'jxc:,— 


,a(x',x'^ 


dX 


/ / / 

■Psjii,AtX^, X^ 


( 23 ) 


Note, that Eq. [22] is different from Eq. (3) of Ref. Jj, where the delta function for the 
solute degrees of freedom is omitted (by mistake?). 


3 Energy representation (ER) 

3.1 Basic definitions in ER 

Collective coordinate. We define a new collective coordinate which is the interaction energy 
between a solute molcule and a solvent molecule: e. We make this coordinate A-independent, 
such that this coordinate is calculated with the solute-solvent potential at full coupling Vsw{'^s,^w), 
irrespective of the ensemble and Hamiltonian which were used to generate this configuration: 

^sji)(Xs,X^) = A=1 (Xg, x^) (24) 


Microscopic density. For a single conhguration of the system the microscopic density in 
energy representation can be written as: 

PL(e) = 5^ ^ (^^««;(x., X^) - e) = (25) 

i=l 


/ + 00 

dx.'Jx.'J{vsU^'s, - e) '^(x', - x,)A(x(„ - x^, 

i=i 




dx^dx^5(v^^{Xs,Xy,) - >^)psw(x^,x^) 


(26) 

(27) 


Potential in ER. The potential in energy representation can be written as: 

/ -l-CX) 

dx^dx^h (n^^(xs, x^) - e) Usnj,x{^s, x^) ( 28 ) 

•CO 

It is important for the following derivation that we choose the lambda path in such a way 
that Usw,x{'^s,^w) is constant on each equi-energy surface of Vsw(x^s,^w)- This can be achieved, 
for instance, when Usw,x{'^s,^w) = Xvswip^sy^w)- With this restriction of the Usw,x potential we 
can write the following identity: 


0 + 00 


i‘sy;p(Xs,Xw) ^ I dfL ■ 6 (VsU^s^X^) - ^) 


(29) 


Taking the partial derivative of the both sides of equation we obtain the formula, which will 
be used later on: 


9Usw,x(X^s^ X^) 

dA 


r*+oo 


de • 6 {vswi^s 


-e) 


9ut^,xi 

dX 


(30) 
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Solute-solvent density distribution in ER. Solute-solvent density distribution in NPT 
ensemble is written as: 



(P(e)). 


Using the definition of microscopic density in ER (Eq. fI7}i and writing explicitly its ensemble 
average we get: 




ir ^ (^) ® Iv x^) - e)p,^(x;, x^) exp [-/3 [/a(x„ x^-)] 


ir ^ (^) ® fv exp [-^Ua(xs, x^»)] 


(31) 


Change of the integration order: 


e f ^/o e ^ [p,^(x,,x^)]exp -^[/a(x,,x(;(™) 

pL,A(e)= / A (v\ -apv r ^ ^ n - r orr) -^ 

J-oo Jo « e j^dx^(ix"“exp [-/3f7A(x^,x"-)J 

The ratio gives us the definition of the solute-solvent density distribution (see Comment 
after Eq. 


^+oo 


pL.Alf) = / 


(32) 


3.2 Kirkwood charging formula in energy representation 

3.2.1 Kirkwood charging formula via density distribution 

Let us obtain the charging formula in energy representation. We start from coordinate repre¬ 
sentation (Eq. [23]): 

Pex = dX J - — -p,^,a(x„x^) = 


Using Eq. [30] we obtain: 

/ I r+oo 

J - 

We change the integration order: 


dx.'Jx.l 


f + CJO 


de6 (u«^(x(,xj - e) 


du' 


sw.X \ 


dX 


Psui,a(x^, X^ 


We use the relation Eq. |32] to obtain: 


dx'Ax'j(v,Ax'„x'j - () p,„AK,K,) 


/".> ,.A 

f^ex — I I de Psw,x\^) 


'0 


dX 


(33) 


Eq. [33] is the Kirkwood’s charging formula in energy representation. 
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3.2.2 Indirect part of potential of mean force (IPMF) 

We can introduce an auxiliary function which is an analogue of the indirect part of 

potential of mean force in coordinate representation: 


PL,A(e) = PL,A=o(e) • exp [-/5 (nL,A(e) + <«,,A(e))] 
The potential then can be rewritten as: 

PL,A(e) 


uL,x{^) =-ksTln 

Psw,X=0\^) 




(34) 


(35) 


3.2.3 Kirkwood charging formnla via IPMF 

Let us rewrite the Kirkwood’s charging formula (Eq. [33]) via 

"1 /■+°° dul^ ;^(e) 

J\ / sw,^\ J 


Hex = d^ de- 


' 0 J —oo 


dX 


-pL,A(e) 


Change of the integration order: 


J —OO J 0 

Integration by parts for the inner integral: 


Hex = I de I dX—— - 


r+oo 

Hex / de 


-1 


PL,A=i(e)“L,A=i(e) - 


j,^PL,A(e) e 
dX --d- u 


dX 


sw,\\ 


Change of the integration order back. Mind that = n|^(e) = e according to the 

definition (Eq. IMland Eq. I28|l : 


, e . ^ , dpUxi^) e , ^ 

Hex = I o?ep.^,A=i(ce - I dX j de --w.^.aIC 


' 0 J —oo 


dX 


(36) 


Let us denote the last term as a functional of the potential and the solute-solvent density 
distribution: 

/•i p+oo dpl„Je) 

J^[pL,xi^),ulw,xi€)]= dX j de ——<^,A(e) (37) 

The functional can be written via IPMF. Here and after, we use the following simplified 
notations: 

g _ g 

Psw,0 Psw,X=0 
Psw Psio,A=l 

Similar notations are adopted for other functions. 

Using Eq. [35] and changing the integration order we obtain from Eq. [371 


-7'[pL,A(e),«L,A(e)] = 


r*+oo pi 

de / dX 

-OO J 0 


dp[ 


sw,X v 


dX 


-ksTln 


Psw,x{ 


Pswfli 


W 


sw,X \ 


( 38 ) 


The first integral in Eq. [38] can be taken analytically by parts: 


dX 


dpLA^) , pL,A(e) 


Ptw,x{^) 


^w,A\ / 1 rsw,A\ / p / \ i rsw.AX / 

Psw,0\^) Psw,0\^) 


/■^ ,,PL,A(e) ^PL,A(e) 

Jo pL,a(^) 
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(39) 


= (U,) In 4^. - {(U^) - pL,o(‘)) 


Psw,o(^) 


as: 


Therefore, we rewrite Eq. [38] using Eq. 

-^[pL,A(e),wL,A(e)] = 


r’+oo 


de 


PsWjOi 


-ksT ( pl„(c) In + f dA-®'’””’"*'' 


'0 


d\ 


-w 


sw,\y 


=)) 


Regrouping the terms we get: 


“+00 


J^[pLAe),ul^.{e)] = kBT / de 


Ptwi^) 


,dpt^Ae) 


(pL(e) - pL,o(^)) - PL(e) In - /5 / 


Psw,o(^) 


dX 


(40) 

This expression can be further simplihed if we choose the A-dependence of the potential 
such that the density distribution is a linear function of A: 


pIw,\{^) — ^ptwi^) + (1 ~ -^)pL,o(f) 
With this restriction (Eq. ITT]) the A-derivative is: 


(41) 


dp[ 


S'W,X ' 


dX 

and the functional (Eq. SOD becomes: 


= (pL(e) -pL,o(e)) 




/ + 00 

de 

•CO 


(pL(e) - pL,o(e)) - PL(e) In - /5(pL(e) - PL,o(e)) [ dXw^,^^^{e) 

Psw,o{^) Jo 

(42) 

Finally, the excess chemical potential (Eq. I5T]) reads: 


"+00 


fe|pLA(e).<.,A(e)l = / ^P^We-UpLaW.^LaW] 


(43) 


3.2.4 Density functional 

For further derivations we would like to consider the functional as a unique functional of 
Psw,xi^)- This can be the case if the solute-solvent interaction potential is a unique functional 
of This implies that there should be only one to which a given Pg^xi^) 

corresponds. Both in coordinate and energy representation it is not the case if we consider 
ensembles where the number of particles is hxed 0,0 [HQ- This can be easily seen from 
the dehnition of ply^^xi^) (Eq- ED and Eq. [2|): if one adds a constant to the solute-solvent 
interaction potential Usw,x the resulting p^^ function does not change (mind, that there is 
a one-to-one correspondence between the potential in energy and coordinate representations 
(Eqs. I2^ and l28|) ). The lack of the one-to-one correspondence between p and u results to the 
fact that the density-density correlation matrix is not invertible and has a singular eigenvalue 

00 [ 3 . 
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Matubayasi proposed a way how to retain the one-to-one p - u correspondence by introducing 
additional condition based on the physical sense. Firstly, he showed that the potentials giving 
different density prohles can differ from each other only by an additive constant (see Appendix 
of Ref. [l^ and Ref. i), Secondly, he set the additive constant to ensure that the chemical 
potential is an intensive property of the system. This can be achieved by ensuring that the 
solute-solvent pair potential reaches zero when particle separation tends to infinity. 

With this approach a one-to-one correspondence between Usw,x{,^) and Psw,\{€.) achieved 
both in coordinate and energy representation. This allows us to consider the potential ;^(e) 
as a functional of Pg^ in a fixed-N ensemble and use the functional calculus to obtain 
approximate free energy functionals. 

Therefore, the excess chemical potential (Eq. I4^ can be written as a density functional of 
the solute-solvent density distribution: 

/ + 00 

dept^{e)e - T[pl^^^{e)] (44) 

■OO 

3.3 Approximate free energy functional. 

The exact free energy functional (Eq. 142]) contains the term which depends on A. To eliminate 
the A-dependence we apply the Percus’s method of functional expansion to obtain approximate 
functionals. 


3.3.1 Hypernetted chain (HNC) - like approximation. 

Following Percus [l^ we obtain the HNC-like approximation by expanding the following func¬ 
tional in powers of density fluctuations pl^id) — p^^ o(^0- 


“+00 




dptw{d) Ptn,i<^')=Ptn, o(^') 

(45) 

With the help of Eq. [35] we rewrite the left hand side of Eq. |35]via IPMF. Therefore, IPMF 
in HNC-like approximation can be written as: 


w 


e,HNC 

sw 


(^) 



de' ■ {p\ 


SW \ 


Psw,0^ 


0) 


d(e-€') dulje) 

PL,o(e) dpiJe') 


Pli^d')=Psn,,o(<^'). 


(46) 


Let us show that the functional derivative in Eq. [46] is a functional inverse of the density- 
density correlation function. For that we start from the dehnition of solute-solvent distribution 
function at full solute coupling: 

pL(e) = (p(e))A=i = 


/o°° d {^) Jy 


(47) 


Let us denote nominator of Eq. |47| as / and denominator as g. Then, End the functional 
derivative of distribution function with respect to solute-solvent potential: 


e e P \ ^9 t 

^PswK^J _ _ 5u1^ _ ^ 

dWgypie") g g g 


(48) 
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Both in / and g only potential energy U depends on To write its explicit dependence 
on we use the relation between the solute-solvent interaction potentials in coordinate and 
energy representations (Eq. [21]); 

17(X^,X^“) = T(Xs) + = 

2=1 

^ + 00 

= T(xJ + ^w,i) - e")^iL(e") + (49) 

i=i 

Therefore, we find the following derivative which will be used in later derivations: 

X rg-/3f/l r+oo 

. e r f\ = -I3e~^^ ^ / de" ■ x^,*) - e")h(e" - e') = 


= '^S{vsnj{^s, ^w,i) - e') = -/3e~^^ pl^{e') (50) 

2=1 

where we used Eq. [281 

With this relation (Eq. [5U]) the first term in Eq. [TS] then can be written as: 

Sf 

^ (51) 

where denotes the ensemble average with the Hamiltonian where the solute-solvent in¬ 

teraction potential is Usw 

Also, with relation Eq. [HI we see that g' = —jdf. With this Eq. [SDlis written as: 

= -13 («„(£'))„„] (52) 

Which equivalently can be written as: 




-P [pL«;(e) e') + pUe)S{e - e') - pL(e)pL(e')] = -/^xL«,(e, e') 


(53) 


where xlww{^i ^') is the density-density correlation function, and Ps^^(e, d) is the three molecule 
distribution density distribution defined by analogy to the two molecule density distribution in 
coordinate representation (see Eq. (2.5.13) of Ref. j^) as: 


P 


e 

sww 


6, e 


/ 

EE 6{v{^,, x^,i) - e)h(n(xs, x^^) 

\ i=i jyi 


From Eqs. [52] and [53] we obtain: 



<^PL(e') 


( ^PL(e) ^ ^ 

UwL(e')/ 


-kBTixLJ-^ (6,e') 


(54) 


(55) 


where (xtww) ^ is fhe functional inverse of the density-density correlation function defined as 
(see Eq. (3.5.8) of Ref. (sf): 


r’+oo 


de" 


■Xs 


(x' 

\A, 


) ^ 
swwJ 




= 6{e 


(56) 
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With Eq. [55] we can rewrite the HNC-like approximation of the indirect part of potential of 
mean force (Eq. H6l) as: 




ksT 


Psw{^) 


PswX) 


(e) 


^+oo 


Psw,Q\ 


*' ■ [pL(e') - PL,ci(e')] ■ (x'L„,ci) ‘ (£.«') 


(57) 


Again, following Percus [13| we obtain the Percus-Yevick-like (PY-like) approximation by ex¬ 
panding the following functional: 


3.3.2 Percus-Yevick (PY) - like approximation. 


“+00 


~ PL,o(e) + / de' ■ (pL(e') “ PL,o(e')) 






( 58 ) 


Psw Psw,Q 


We rewrite Eq. EHl via IPMF (Eq. 


w 


e,PY 


[e) = —/cfiTln (l + J de' ■ 




PL,o(e) 


Psw Psw,Q^ 

With the help of Eq. |36]we can rewrite the PY-like approximation via the HNC-like w. 

= -ksTln (1 - /3wf^^^{e)) ( 60 ) 

3.3.3 Lambda-integral in HNC-like approximation. 

When is the solute-solvent interaction potential the corresponding IPMF is written as: 


(59) 


e.HNC/ \ 1 rji 

Wsw,x (e) = -kpT 


Plw,x{P) - pL,o(^) 
PL,o(e) 


^H-oo 


*' ■ [ptw,\W) - /’L,o(e')] ■ {xLwp ' (£. e') 


(61) 

With the linear dependence of x E]) Eq. ED can be written via at full 

solute coupling: 

^ (62) 

The A-integral in Eq. [69] can be written in HNC-like approximation as: 


P I = wfr^ie) ■^j^dX-X = 

3.3.4 Lambda-integral in PY-like approximation. 


( 63 ) 


When Ux is the solute-solvent interaction potential the corresponding IPMF in PY-like approx¬ 
imation is written as (see Eq. [6U]) : 

(e) = -kpT In (l - ldwl'^x^{e)^ = -keT In (l - A • (64) 

Now, we use the following known tabulated relation: 

, , , ,, (ax + b)-\n(ax + b) — ax 

/ dx ■ ln(ax + b) = - 
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to find the A-integral: 


/ 3 / dXw:^Z{e) = 


+ l] • In ^^(e) + l] + /3wfjf^^{e) 


(65) 


Next, we use the relation between w in PY and HNC-like approximations at full solute 
coupling A = 1 (see Eq. | 


= -ksTin (1 - => = e-'’”™ '<> - 1 (66) 

Using Eq. [66] we rewrite Eq. 


as: 


Id / dXwl’^^ie) = - In [1 - /3<;f^(e)] + 1 + 


In [l - I3wfj^^{e)] 


(3w 


e,PYf N 
sw i ^ j 


(67) 


3.3.5 Constructing hybrid functional. 

The approximate functional is developed by Matubayasi and Nakahara based on the following 
considerations. To make an end-point expression of fiex we need to approximate the A-integral 
in Eq. 02] Beforehand, we would like to note that the approximate expression of the A-integral 
combines four parts. 

Firstly, the A-integration can be analytically performed both in PY-like and HNC-like ap¬ 
proximations (see Eqs. 1671 and l63|l . There is a general knowledge in the field that in the case 
of simple liquids the PY approximation works better for short range repulsive potentials, while 
HNC approximation performs better for long-range attractive potentials j8|. Matubayasi and 
Nakahara decided to use the PY-like expression for the A-integral in the unfavorable region 
of solvation (tc®^ > 0) and HNC-like expression for the A-integral in the favorable region of 
solvation (tc®^ < 0). 

Secondly, the indirect part of potential of mean force can be determined from unbi¬ 
ased molecular simulations (regular MD or Monte-Carlo) only outside of the solute-core region 
(region of very large solute-solvent interaction energies: e). Matubayasi and Nakahara pro¬ 
posed to use HNC-approximation of the potential of mean force when is not resolved from 
molecular simulations. In HNC-like approximation is determined by the solute-solvent 

density distribution pl^ q and the inverse of the density-density correlation function {xtww o) 
in the case of the zero solute-solvent coupling (this can bee seen from Eq. (57] where the differ¬ 
ence — Ps^fi can be safely approximated by —o since in the core region pl^ « pI^^q). 
Therefore, in the core-region can be calculated with high resolution in the ensemble, 

where solute and solvent are fully decoupled and the probability to hnd solvent molecule in the 
core region is high. The later can be in the most convenient way realized by the insertion of the 
solute molecule conhgurations into the ensemble of precalculated pure solvent conhgurations. 

Combination of the two different expression for the A-integral and the two different input 
w functions results into functional consisting of four parts. The final expression for the excess 
chemical potential (Eq. HT|) is: 


“+00 


/^ex[pL(e),pL.o(e)>xL«,,o(ce')] = / “-^[pL(e)> dL,o(e), xL«,,o(c e')] (68) 


where 
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“+00 


= kBT 


de 


(pL(e) -pL,o(e)) -pL(e) In 


PL(e) 


- (pL(e) -pL,o(e)) •2:[p sw) Psw^O") 

( 69 ) 


where X is the approximated A-integrah 


/5 / X[p^^, xLwfl] = «(e) • Fw{e) + [1 - a(e)] • F^HNc{e) (70) 

Jo 

where the functions F,„ and F,„hnc are in turn written as the combination of PY and HNC- 


like expressions for the A-integral: 

I PwlwiF) 

FUe) = ' 


2 


— 1' 


when > 0 

when wl^{e) < 0 


(71) 


and 




^ when > 0 

FwHNc{Fj = { ^ In (l - Bw'^FNCf (72) 

' L + ^ ^ when w^^NCr a ^ q 

^ „ e,HNCf \ ’ wnen ^ u 

pWsw (e) 

The parameter a(e\ which is responsible for merging parts with different w functions, is 


- In [1 - /3<^^^^(e)] 

ter Q:(eh ^ 
set heuristically as jll: 


a(e) = 



PL(e) -PL,o(e) 
Psw(^) + Ptw,o(^) 


when p^^(e) > pL,o(e) 
, when p^^(e) < pL,o(e) 


(73) 
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